【直播】我的基因组20：覆盖度详细探究
前面我们在第8讲提到了公司给我的一个报告的统计表格，有人反映不会做。
本来应该只需要给我6亿条reads的（PE150测序，人30X），但是足足给了我8.9亿条！（但事实上很多paper发表的基因组高于60X的也不少）
表格里面提到了好几个概念，比如duplicate的reads，一般说的PCR造成的duplicate，在找变异的时候需要去除掉。然后是那些比对到了不同染色体的reads pair，虽然只有2.29% ，也是需要重点分析的。（前面我也讲了如何提取以及分别分析它们！）
如果只是想重新前面的这些统计指标，非常很简单，就是samtools工具提供了一个flagstat功能，用法及结果如下：
samtools flagstat P_jmzeng.final.bam
899361748 + 0 in total (QC-passed reads + QC-failed reads)
8597742 + 0 secondary
0 + 0 supplementary
132556557 + 0 duplicates
890858540 + 0 mapped (99.05% : N/A)
890764006 + 0 paired in sequencing
445382003 + 0 read1
445382003 + 0 read2
853255862 + 0 properly paired (95.79% : N/A)
881249604 + 0 with itself and mate mapped
1011194 + 0 singletons (0.11% : N/A)
20382234 + 0 with mate mapped to a different chr
12511988 + 0 with mate mapped to a different chr (mapQ>=5)
从结果很明显可以看到，公司的确给了我8.9亿的reads，这个没错，duplicate的情况也的确是132556557/899361748 = 14.73%，而mapping的情况和proper mapping的情况也都显示好了，可以看到公司用的是同样的命令和方法！没有什么神秘的，我们生信工程师，做得就是这个，而且可以做得更好。
后面的探索全基因组区域中碱基覆盖深度不低于多少X的比例，是需要画一个图，有非常多的现成的工具可以使用，包括 、 ，还有的几个命令。
其实本身原理很简单，就是把全基因组的每个坐标的depth都得到，然后得到depth的频数，然后画图。
我们可以对每条染色体单独来绘图，也可以针对全基因组来绘图，当然，公司给我们的数据是针对全基因组的。 脚本很简单：
samtools mpileup P_jmzeng.final.bam |perl -alne '{if($F[3]>100){$depth{"over100"}++}else{$depth{$F[3]}++}}END{print "$_\t$depth{$_}" foreach sort{$a <=> $b}keys %depth}' 当是需要运行很久，毕竟全基因组的bam文件太大了。
如果分开运行，可以对下面的各条染色体bam文件批量跑一个脚本
脚本如下：
ls P_jmzeng.final.REF*.bam |while read id
do
echo $id
samtools mpileup $id |perl -alne '{if($F[3]>100){$depth{"over100"}++}else{$depth{$F[3]}++}}END{print "$_\t$depth{$_}" foreach sort{$a <=> $b}keys %depth}' >$id.depth.txt
done
跑完之后，对每条染色体都会输出如下文件：
over100 110789
0 27065
1 1730286
2 2219409
3 2728526
4 3251046
5 3774335
6 4303971
~~~~~~~~~~~~~~~~~~~~~~~
后面省略94行，每一行都是两列，第一列是测序深度，第二列是有着该测序深度的位点是多少个！所以行的第二列加起来，就是染色体的长度！！！面这个例子是X染色体的，绘图如下，可以看到X染色体的测序深度其实并不怎么好，全基因组测序深度平均高达44X，可是这个X染色体超过44X的只占极少的比例。
我用excel表格简单画个图如下：（当然，作为一个高级生物信息学工程师，用excel表格是有点low，但是这里只是为了说明一个问题，我们后面还是写程序的，用R语言）
比如我们再看看10号染色体，我随意在R里面画了个图：
a=read.table('P_jmzeng.final.REF_10.bam.depth.txt',stringsAsFactors = F)
plot(a[,2],type = 'l',xaxt="n",lwd=3)
axis(1, at = 1:nrow(a),labels =a[,1] , las=1)
abline(v=44)
可以看到10号染色体的测序深度要显著好于X染色体，大部分的测序深度都超过了30X！！（不过，这里虽然用来R，但是出图很丑，虽然不是专业做可视化的，但是想调整美观一点问题也不大，就是好耗费时间。）
做出depth的累积曲线图了，也是很简单的！
a=read.table('P_jmzeng.final.REF_10.bam.depth.txt',stringsAsFactors = F)
over100=a[1,]
a=a[-1,]
a=rbind(a,over100)
a[,3]= cumsum(a[,2])
a[,4]= 1-a[,3]/sum(a[,2])
plot(a[,4],type = 'l',xaxt="n",lwd=3)
axis(1, at = 1:nrow(a),labels =a[,1] , las=1)
abline(v=10);abline(v=20);abline(v=30)
最后得到的数据如下：第4列，就是不低于多少X的比例
可以看到大于10X测序深度的比例仍然高达92.88%，效果杠杠的！！
这个测序深度累积分布图，就是很多人的重点！！！！（请大家仔细学习我上面的统计脚本和画图方法）
加油吧，骚年！！！！
参考：
文：Jimmy、阿尔的太阳
图文编辑：吃瓜群众